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Abstract 

We investigate two-dimensional Ising systems with multisping in- 
teractions of three- (m = 3) and four-body terms (m = 4). The 
application of a new type of finite-size algorithm of de Oliveira allow 
us to clearly distinguish a first-order transition (in the m = 4 case) 
from a continuous one (in the m = 3 one). We also study the damage 
spreading in these systems. In this study, a dynamical phenomenon is 
observed to occur at a critical point separating a chaotic phase from 
a frozen one. However, the width of the interval where this transition 
happens does not yield a conclusive evidence about the order of the 
phase transition. 
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1 Introduction 



The presence of multispin interactions can lead to severe alterations in the 
behavior of statistical mechanical models such as the Baxter and the Ashkin- 
Teller models. There are other models of the same type and these can include 
three-, four- or more-body terms in dimension d > 2. For instance, we may 
cite the Ising system with multispin interactions. For a two-dimensional 
problem, the Hamiltonian H of such an Ising-like system is expressed by 

1 ( m—l ~\ 



l=0 



where ks is the Boltzmann constant, T is the temperature and Sij = ±1 is 
a boolean variable localized on site 

Turban and Debierre Jl| have looked at this anisotropic model and shown 
that it has a single transition located at the self-dual critical point, once 
the above Hamiltonian is self-dual for any m. The critical point can be 
determined through the known relation sinh(2i^ :r ) sinh^fCy) = 1, which is 
independent of m. For the isotropic case [K x = K y ), the physical solution 
of this relation is K x = K y = K c = |ln(l + V2) = 0.44068679... . From 
the spin reversal invariance of the Hamiltonian as well as from an analysis of 
the energy of domain walls, one can see that the ground state is 2 m_1 times 
degenerate and expect that the model will be on the same universality class 
as the g-state Potts model whenever q = 2 m_1 [1,2]. 

Mean-field theory [[[]] and finite size scaling [1,3-6] were first used to study 
this problem. Further improvements of the mean-field methods and Monte 
Carlo simulations |J helped to indicate a first-order transition in the m = 4 
case. The method of the fourth-order cumulant of Challa, Landau and Binder 
and a Monte Carlo histogram technique of Ferrenberg and Swendsen [T^] 



were then used to study the order of the transitions in the m = 3 and m = 4 
cases 0: the phase transition was characterized to be continuous in the 
m = 3 case and was asserted to be a first-order one in the m = 4 case. 

In the present study, we apply a new type of finite-size scaling algorithm 



TTf to study the critical properties of the above model. We clearly distinguish 
a first-order transition (in the m = 4 case) from a continuous one (in the 
m = 3 case). In addition, we obtain good collapsed curves when the pertinent 
exponents are used. 

We also consider the dynamics of spreading phenomena [12,13] in the 
system. A critical parameter which separates a chaotic phase from a 
frozen one is obtained in both cases (m — 3 and m — 4). However, unlike 
the finite-size scaling method, a damage spreading analysis [14-16] seems to 



2 



be not able to characterize the order of the phase transition. 



2 Finite-Size Scaling 

The finite-size scaling algorithm of de Oliveira [11,17] is based on two thermo- 
dynamic quantities, namely the bulk quantity Q and the surface correlation 
function r. 

For an Ising system, r can be calculated as follows: one considers two 
opposite surfaces of the lattice and verifies which spin state dominates each 
one of them. A counter is increased by 1 if both surfaces are dominated by 
the same state, otherwise it is decreased by 1. This counter is normalized 
by the number of measurements and the result is denoted by r. In this way, 
r is a step function of T — T c . For T < T c the surfaces are well-correlated 
and r = 1, whereas for T > T c the surfaces are uncorrelated and thus r = 0. 
This criterion holds for second-order transitions where the correlation length 
diverges; however, at a first-order transition the correlation length remains 
finite and thus a multidomain state is possible even below T c , which implies 
that this normalized counter may be close to zero in some samples and close 
to 1 in other samples fllTf . 

The bulk quantity Q is the average of the sign of the sum of the Ising 
spins ±1 [11,18]. 

Both quantities Q and r are shown [11,18] to scale as L° at the criti- 
cal temperature T c and in the thermodynamic limit (L is the linear size of 



the system). So they behave like Binder's fourth-order cumulant [19 .The 



function Q is based on a bulk measure of the majority of spins, whereas r 



measures the correlation between two opposite surfaces of the system [17| 



Another important quantity is the bulk magnetization M. For an Ising 



system, it is the average of the absolute value of the sum of all spins |TT 
At the critical temperature, this quantity scales as M ~ L y , where y is 
the magnetic exponent. For a continuous phase transition, the exponent 
is y = D — where D is the geometrical dimension of the system, j3 is 
the exponent associated to the spontaneous magnetization (order parameter) 
and v is the thermal critical exponent which governs the divergence of the 
correlation length. The definition of y can be extended even to a first-order 
transition, where the exponent v cannot be defined. In this case, M ~ L D 
and thus y = D 0. 

We simulate two-dimensional Ising systems with multispin interactions 
and Hamiltonian expressed by (|]) for the cases m = 3 and m — 4. For each 
Monte Carlo updating, spins are selected sequentially and flipped with the 
thermal probability of the Metropolis algorithm. Simulations are performed 
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on the square lattice with periodic boundary conditions for different lattice 
sizes L. For fixed K x (= K y ), we average the quantities r and M over a total 
of (2 — 4) x 10 5 Monte Carlo steps after discarding the initial 4 x 10 4 transient 
steps. In order to obtain these averages, we consider a subset of spins be- 
longing to the rows of sites which are separated from each other by multiples 
of m lattice parameters. In this way, we may calculate the magnetization 
as the average of the modulus of m times the sum of all spins belonging to 
those rows. The surface correlation function may be calculated as defined 
above provided that the two mentioned opposite surfaces correspond to two 
of those rows which are L/2 lattice parameters apart from each other. The 
statistics is improved by averaging over all possible subsets of rows of sites. 

The bulk magnetization M is measured at K c for different linear sizes L. 
The plots of M versus L in logarithmic scale are presented in Figure 1 (for 
simulations with m — 3) and in Figure 2 (for m = 4). The straight lines 
confirm the scaling relation M ~ L y . For m = 3, the magnetic exponent 
obtained is y = 1.83 ±0.03. This result is really consistent with the expected 
second-order value y = D — - = 2 — A = 1.8125 (this model is considered 



to be on the same universality class as the q = 4 Potts and Baxter- Wu |20 
models with exponents a = v = 2/3 but (3 = 1/8 [0]). For m = 4, we 
obtain the value y = 2.08 ± 0.05 which points to the expected first-order 
transition with magnetic exponent y = D = 2 (remind that the m = 4 case 
should be on the same universality class as the q = 8 Potts model). So this 
finite-size scaling method seems to be still very trustworthy when applied 
to Ising systems with multispin interactions. It is a way to provide a clear 
determination of a first-order phase transition as well as a second-order one. 



Figure 1 to be inserted here. 



Figure 2 to be inserted here. 



Additionally, plots of the surface correlation function r versus L 1 ^ (K x — 
K c ) for different sizes L can lead to estimatives of the critical coupling K c 
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and the thermal exponent v by adjusting these parameters so that all points 
collapse onto the same curve |I7| . We instead prefer to show that we do 
obtain this collapse if the known values of K c and v are previously used. 
Data from simulations with m = 3 lead to the plot of r versus L 1 ^ (K x — K c ) 
(with v = 2/3) as shown in Figure 3. In this graph, the points correspond 
to discrete values of K x nearby K c for lattice sizes L = 30 (+), 60 (*) and 
90 (□). A good collapse is observed for all three sets of points. In the case 
m = 4 where the thermal exponent v is not defined we can also obtain a 
collapse if we consider the artificial value v — 1/2. In Figure 4, we thus 
plot t versus L 2 (K x — K c ) for m = 4 and lattice sizes L = 32 (+), 64 (*) 
and 96 (□). One should consider the standard deviations At of some points 
when appreciating the accuracy of the collapse: in our simulations statistical 
errors may reach the value Ar ~ 0.1 for those points on the fall of curve 
corresponding to L = 96 (□). 



Figure 3 to be inserted here. 



Figure 4 to be inserted here. 



3 Damage Spreading 

Of late years the damage spreading method has been used as a numerical 
approach to study the propagation of a perturbation throughout many sys- 
tems such as spin glasses [fT2"| , Ising model [13-15], g-state Potts model [16[ 



and cellular automata In the present work, we study damage spreading 
in Ising systems with multispin interactions. We again investigate the Ising- 
like system (on the square lattice) with a Hamiltonian described by equation 
(JJ) with K y = K x . Consider that the time evolution of two independent spin 
configurations of the system is governed by the same dynamics and the same 
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sequence of random numbers generated in the Monte Carlo process. Let 
A t = {S£j(t)} and B t = {S^(t)} be these two configurations at time t. The 
total damage D(t) between them is defined as the fraction of corresponding 
spins with different signs, that is, 

D ^ = ^ Ei^lw-sgwi, ( 2 ) 

where N is the number of lattice sites. 

A starting configuration is thermalized at a fixed K x (oc l/T) by the 
Glauber dynamics. At time t — 0, this thermalized configuration is termed 
Aq and a second configuration Bq is created from Aq by flipping ("damag- 
ing") a chosen fraction D(0) = M/N of corresponding spins. Then both 
configurations evolve in time according to this same dynamics and with the 
same sequence of random numbers. After a transient (2000 steps per site), 
we get a time average of D(t) (over 7000 Monte Carlo steps). We calculate 
its medium value < D(t) > over several samples of configurations A and 
B , for each K x and initial damage D(0). From a total of 50 initial samples, 
we only average over those configurations where damage is non null (unless 
all of them present a damage that has become equal to zero). 

In a Monte Carlo step at time t, all lattice sites are sequentially 

visited and each spin Sij(t) is flipped with probability 

where AH is energy change associated with such a possible spin-flip. The nu- 
merical procedure for updating the spins consists in generating a sequence of 
random numbers rij(t) uniformly distributed in the interval [0, 1] and setting 
Sij(t+1) = —Sij(t) iirijit) < Pij(t) or setting Si j (£+1) = Sij(t) otherwise. 
One must use the same sequence r iy j{t) for updating both configurations A t 
and B t . 

We simulate Ising-like systems on square lattices L x L with periodic 
boundary conditions. The averaged long-time damage < D(t) > versus K x 
for the case m = 3 (with L = 42) is presented in Figure 5 whereas results 
corresponding to m = 4 (with L = 40) are shown in Figure 6. For each case 
we plot two curves: one of them connects the points (squares) obtained from 
simulations with initial damage D(0) = l/N = 1/L 2 and the other one (star 
points) corresponds to D(0) = 0.90. 
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Figure 5 to be inserted here. 



Figure 6 to be inserted here. 



In both graphs, for the lower values of K x (corresponding to the higher 
temperatures) the long-time damage reaches the value D* = 1/2 which is the 
same result from Glauber dynamics for damage spreading in the Ising Model. 
In fact, there are also two states per spin in the multispin Ising system so 
two corresponding spins S^At) and S?At) present 4 possible configurations. 
For those values of K x , all 4 possibilities should be equally probable but 
only half of them would correspond to damaged configurations. So there is a 
range of low values of K x where the damage does not depend on the initial 
conditions. However, exactly at K x = 0, one should note Pi,j(t) = 1/2 so 
that two corresponding spins would (or would not) simultaneously change 
their states thus preserving the initial damage (< D(t) >= D(0)). There is 
also another situation where the initial damage is preserved: for D(0) = 1 
and m = 4 one can easily prove that two corresponding spins have the same 
probability Pij(t) of flipping (the preservation of initial damage D(0) = 1 
has been already observed for the m = 2 Ising ferromagnet). 

For higher values of K x (or lower temperatures) the damage depends on 
the initial conditions. In this case, we observe that if the initial damage is 
D(0) < 1/2 then the long-time damage is < D(t) >= for both m = 3 and 
m = 4. On the other hand, if D(0) > 1/2 then < D(t) >^ 2/3 for m = 3 
whereas < D(t) >— > 1 for m = 4. This behaviour is illustrated in Figures 5 
and 6 for the chosen values D(0) = 1/L 2 < 1/2 and D(0) = 0.90 > 1/2. 

For increasing K x close to K c the damage changes from D* — 1/2 to the 
value corresponding to the second plateau of each curve. Oscillations in the 
critical region are due to statistical fluctuations. The standard deviations of 
those points in the critical interval may vary from ~ 0.05 unit (star points) 
up to ~ 0.15 (squares). 

Similar dynamical critical phenomenon has been already observed in the 
literature for Ising jnj and g-state Potts models (with Hamiltonian given by 
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H = —JYl<ij>fi( 'ii (7 j)) EH- m reference |IE| , it has been shown that sim- 



ulations of damage spreading in the g-state Potts model through a Glauber 
dynamics yielded a frozen phase for a temperature T < Td and a chaotic 
phase for T > T^. Damage was observed to assume a high-temperature 
value {q — l)/q above a characteristic temperature T* > T^. This allowed the 
authors to define an interval AT = T* — Td whose q dependence presented a 
decreasing behaviour with increasing values of q, giving "some indication of 



the order of the phase transition" fll6 



The above fact has motivated us to investigate the order of the phase 
transiton in the m = 3 and m = 4 multispin Ising systems through a damage 
spreading analysis. Our results above (expressed as functions of the coupling 
K x instead of T) show that, in fact, there is a dynamical critical phenomenon 
at a critical coupling Kd so that for K x < Kd the damage is insensitive to 
its initial value {chaotic phase) whereas for K x > there are two definite 
constant values of damage corresponding to initial conditions D(0) < 1/2 and 
D(0) > 1/2 (frozen phase). We estimate Kd = 0.440(1) for m = 3 (Figure 
5) and Kd = 0.444(1) for m = 4 (Figure 6). We also identify a characteristic 
coupling K* < Kd so that for K x < K* the damage assumes the value 
D* = 1/2 and the curves corresponding to different initial conditions join 
themselves. Then we estimate an interval AK = Kd — K* = 0.006(2) for m = 
3 and AK = 0.008(2) for m = 4. The first result might be compared with 

the interval AK ~ ^ (± - ±) ~ |(ak ~ ok) ~ °- 006 ( for the g = 4 
Potts model) extracted from Fig. 3 of reference |16[ , where the non vanishing 
interval was taken as an indicator of the continuous phase transition. On the 
other hand, the estimative of the interval AK for m = 4 does not allow us 
to characterize the order of the phase transition (which is discontinuous in 
this case). 



4 Conclusion 

In summary, in order to investigate the phase transition of Ising systems 
with multispin interactions we have applied two methods: a finite-size scaling 



algorithm of de Oliveira [£L1[] and a damage spreading analysis. 

Through the finite-size scaling method we have obtained the scaling re- 
lation for bulk magnetization, M ~ L y , with magnetic exponents consistent 
with the expected second-order transition value y = 1.8125 (for m = 3) and 
first-order transition value y = D = 2 (for m = 4). In this way the orders of 
phase transitions in the m = 3 and m = 4 multispin Ising systems have been 
clearly distinguished. Additionally, plots of the surface correlation function 
t versus L X I V [K X — K c ) for different sizes L have collapsed onto the same 
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curve, showing that this technique can also be used as a way to estimate K c 
and v. 

Regarding the analysis of damage spreading in the system, we have ob- 
served a dynamical critical phenomenon at a critical coupling m K c sepa- 
rating a chaotic phase (where damage is insensitive to its initial value) from 
a frozen phase (with two definite constant values of damage). However, the 
width AK of the interval where this transition occurs has not yielded a 
conclusive evidence about the order of the phase transition. 
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FIGURES CAPTION 



Figure 1 - Bulk magnetization M at K c from simulations with m = 3 for 
linear lattice sizes L = 30, 42, 54, 66, 78, 90 , 120 and 150. The straight line 
confirms the scaling relation M ~ L y with magnetic exponent y = 1.83±0.03 
which is consistent with the value y = D — £ = 1.8125 expected for a second- 
order transition. 

Figure 2 - Bulk magnetization M at K c from simulations with m = 4 for 
linear lattice sizes L = 32, 48, 64, 80, 96, 112 , 128 and 144. The straight line 
confirms the scaling relation M ~ L y with magnetic exponent y = 2.08±0.05 
pointing to the expected first-order transition value y = D = 2. 

Figure 3 - Plot of the surface correlation function r versus L 1 ^ (K x — K c ) 
(with v = 2/3) for m = 3 and lattice sizes L = 30 (+), 60 (*) and 90 (□). 

Figure 4 - Plot of the surface correlation function r versus L 2 (K x — K c ) 
for m = 4 and lattice sizes L = 32 (+), 64 (★) and 96 (□). 

Figure 5 - Average damage < D(t) > versus K x for the m = 3 multispin 
Ising system with L = 42. The squares (□) represent points obtained from 
simulations with initial damage -D(O) = 1/N = 1/L 2 < 1/2 and the star 
points (*) correspond to -D(O) = 0.90 > 1/2. 

Figure 6 - Average damage < D(t) > versus K x for the m = 4 multispin 
Ising system with L = 40. The squares (□) represent points obtained from 
simulations with initial damage -D(O) = 1/N = 1/L 2 < 1/2 and the star 
points (*) correspond to D(0) = 0.90 > 1/2. 
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FIGURE 1 
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FIGURE 2 
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FIGURE 3 
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FIGURE 4 
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FIGURE 5 
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FIGURE 6 
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